#include "stdafx.h"


#include <stdio.h>

#include "EphemEPM.h"
#include "CustomException.h"
#include "EphemWrapper.h"
#include "SimpleFile.h"
#include "AsteroidOrbitParam.h"
#include "Trajectory.h"
#include "Approach.h"
#include "RadaParamInterface.h"
#include "DEEphem.h"

struct SFileData
{
	double dDate;
	double dVal;
	double dCalcVal1, dCalcVal2;
	int nCenter;
	int nTarg;
	int nPos;
	int nLine;
	int nCheck;

	SFileData()
		: dDate(0.0)
		, nCenter(0)
		, nTarg (0)
		, nPos (0)
		, dVal (0.0)
		, dCalcVal1(0.0), dCalcVal2(0.0)
		, nLine(0)
		, nCheck(0)
	{};

	bool Read(FILE * file, int line)
	{
		int nwasRead = fscanf(file, "%*d %*10c %lf %d %d %d %lf", &dDate, &nTarg, &nCenter, &nPos, &dVal);
		--nTarg;
		--nCenter;
		--nPos;
		nLine = line;
		return (nwasRead == 5);
	}
};

int CheckDEephem()
{
	//CEphemWrapper ephem("../Data/epm2008", CEphemWrapper::eEPM);
	CEphemWrapper ephem("../Data/de423.dos", CEphemWrapper::eDE);

	CDeEphemeride de("../data/de423.dos");
	//	de.GetCoord(eMERCURY, 2378496.5000000000);

	int nLine = 1;
	int nError = 0;
	int nCheck = 0;
	CSimpleFile testFile("testpo.423", "r");
	double dSummError1 = 0.0;
	double dSummError2 = 0.0;


	typedef list<SFileData> ListTestData;
	ListTestData listTestData;

	while(! feof(testFile))
	{
		nLine ++;
		SFileData testData;
		if(! testData.Read(testFile, nLine) )
			break;

		if( testData.nTarg > 10 || testData.nCenter > 10)
			continue;
		listTestData.push_back(testData);
	}

	int t1 = 3, t2 = 3;
	DWORD dwStart = GetTickCount();
	int ntry = t1;

	while(ntry --> 0)
	{
		for(ListTestData::iterator it = listTestData.begin(); it!=listTestData.end(); ++it)
		{	
			nCheck ++;
			SFileData & testData = *it;
			int nTarg = testData.nTarg;
			int nCenter = testData.nCenter;
			int nPos  = testData.nPos ;
			testData.nCheck = nCheck;

			bool bBaryCenter = true;
			bool bGetVelocity = true;

			AllBodies coord;
			int tt = t2;
			while ( tt --> 0)
				ephem.GetCoord(testData.dDate, bBaryCenter, coord, bGetVelocity);

			double dResVal = 0.0;

			if(nPos < 3)
				dResVal = coord.bodies[nTarg].x[nPos] - coord.bodies[nCenter].x[nPos]; 
			else
				dResVal = coord.bodies[nTarg].v[nPos-3] - coord.bodies[nCenter].v[nPos-3]; 

			testData.dCalcVal1 = dResVal;
		}
	}
	dwStart = GetTickCount() - dwStart;

	DWORD dwStart2 = GetTickCount();
	int ntry2 = t1;
	while(ntry2 --> 0)
	{
		for(ListTestData::iterator it = listTestData.begin(); it!=listTestData.end(); ++it)
		{	
			nCheck ++;
			SFileData & testData = *it;
			int nTarg = testData.nTarg;
			int nCenter = testData.nCenter;
			int nPos  = testData.nPos ;
			testData.nCheck = nCheck;

			bool bBaryCenter = true;
			bool bGetVelocity = true;

			double dd[12][6];
			int tt = t2;
			while ( tt --> 0)
				de.GetCoord1(testData.dDate, dd);

			double dResVal2 = .0;

			dResVal2 = dd[nTarg][nPos] - dd[nCenter][nPos]; 

			testData.dCalcVal2 = dResVal2;
		}
	}
	dwStart2 = GetTickCount() - dwStart2;

	printf("===== %d %d =======\n", dwStart, dwStart2);

	for(ListTestData::iterator it = listTestData.begin(); it!=listTestData.end(); ++it)
	{	
		SFileData & testData = *it;
		double dDiff1 = abs(testData.dCalcVal1 - testData.dVal);
		double dDiff2 = abs(testData.dCalcVal2 - testData.dVal);
		double DRelDiff1 = 0.5 * (dDiff1) / (testData.dCalcVal1 + testData.dVal);
		double DRelDiff2 = 0.5 * (dDiff2) / (testData.dCalcVal2 + testData.dVal);
		// 		dSummError1 += DRelDiff * DRelDiff;
		// 		dSummError2 += fabs(DRelDiff);

		if(0)
			printf("%d\t%2d %2d %2d  %+25.20lf %+25.20lf %+20.10e %+20.10e %c\n", 
			testData.nLine, testData.nTarg+1, testData.nCenter+1, testData.nPos+1, testData.dCalcVal1, testData.dVal, dDiff1, DRelDiff1, (testData.nPos < 3)? 'c': 'v');

		double eps = 1.0e-14;
		if(fabs(dDiff1) > eps || fabs(dDiff2) > eps)
		{
			printf("== error at line %d (%d %d %d) \t(diff %20.10e, %20.10e | %20.10e)\n", 
				testData.nLine, testData.nCenter, testData.nTarg, testData.nPos,
				testData.dCalcVal1 - testData.dVal, testData.dCalcVal2 - testData.dVal, testData.dCalcVal1 - testData.dCalcVal2);
			nError ++;

		}
		//printf("\n");
	}

	// 	dSummError1 = sqrt(dSummError1);
	// 	printf("\n\n %d line was read, was %d check, %d error found \nErr sum %e %e (%e %e) %d ms\n", 
	// 		nLine, nCheck, nError, dSummError1, dSummError2, dSummError1/nCheck, dSummError2/nCheck, dwStart);

	printf("Press any key...");

	getc(stdin);
	return 0;
}

int CompareEmpDe()
{
	CEphemWrapper ephemEPM("../Data/epm2008", CEphemWrapper::eEPM);
	CEphemWrapper ephemDE("../Data/de423.dos", CEphemWrapper::eDE);

	double date = 2455400.5;
	bool bBaryCenter = true;
	bool bGetVelocity = true;

	AllBodies coordEPM;
	AllBodies coordDE;

	ephemEPM.GetCoord(date, bBaryCenter, coordEPM, bGetVelocity);
	ephemDE.GetCoord(date, bBaryCenter, coordDE, bGetVelocity);

	for(int bd = 0; bd < BODY_COUNT; ++bd)
	{
		printf("\n%d:  ========================= %s ========================= \n", bd, CEphemWrapper::GetBodyName((eBodies)bd));

		SBody & bodyEPM = coordEPM.bodies[bd];
		SBody & bodyDE = coordDE.bodies[bd];

		printf("GM \t%+20.15lf \t%+20.15lf \t%+20.15e\n\n", bodyDE.GM, bodyEPM.GM, (bodyDE.GM - bodyEPM.GM) * 2.0/(bodyDE.GM + bodyEPM.GM) );

		for(int i = 0; i < 3; ++i)
		{
			printf("x[%d] \t%+25.20lf \t%+25.20lf \t%+10.5e \t%+20.15lf km\n",   i, bodyDE.x[i], bodyEPM.x[i], 
				(bodyDE.x[i] - bodyEPM.x[i]) * 2.0/(bodyDE.x[i] + bodyEPM.x[i]),
				CEphemWrapper::AUtoKM(bodyDE.x[i] - bodyEPM.x[i]));
			printf("v[%d] \t%+25.20lf \t%+25.20lf \t%+10.5e \t%+20.15lf m/s\n\n", i, bodyDE.v[i], bodyEPM.v[i], 
				(bodyDE.v[i] - bodyEPM.v[i]) * 2.0/(bodyDE.v[i] + bodyEPM.v[i]),
				CEphemWrapper::AUDtoMS(bodyDE.v[i] - bodyEPM.v[i]) * 1000.0);
		}
	}

	printf("Press any key...");
	//getc(stdin);
	return 0;
}

int CheckEpm(int argc, char* argv[])
{
	double dj;          // Julian date (frac part)
	int status;         // success
	double xv[3], vv[3];// positions and velocities
	int pl_num;         // number of object
	int centr_num;      // number of center
	FILE *out;

	CEphemEPM epm("../Data/epm2008");

	char out_file[80] = "\0"; // result file

	dj = 2378500.5;

	sprintf(out_file, "%10.2f", dj);
	out = fopen(out_file, "wb");

	fprintf(out, "%s", "The numbering convention for 'PL_NUM' and 'CENTR_NUM':\n");
	fprintf(out, "%s", "                1 = MERCURY           8 = NEPTUNE\n");
	fprintf(out, "%s", "                2 = VENUS             9 = PLUTO\n");
	fprintf(out, "%s", "                3 = EARTH            10 = MOON\n");
	fprintf(out, "%s", "                4 = MARS             11 = SUN\n");            
	fprintf(out, "%s", "                5 = JUPITER          12 = SOLAR-SYSTEM BARYCENTER\n");                
	fprintf(out, "%s", "                6 = SATURN           13 = EARTH-MOON BARYCENTER\n");                
	fprintf(out, "%s", "                7 = URANUS\n");                        

	fprintf(out, "%10.2f\n", dj);
	printf("%10.2f\n", dj);
	for(centr_num = 1; centr_num <=13; centr_num++)
	{
		fprintf(out, "%s %i\n", "Coordinates regard to ", centr_num);
		for(pl_num = 1; pl_num <= 13; pl_num++)
		{
			// To calculate position and velocity of object 'PL_NUM'
			// with respect to the object 'CENTR_NUM' at the date 'jd.dj'
			// call subroutine calc_EPM(pl_num, centr_num, jd, dj, xv, vv, status)
			status = epm.calc_EPM(pl_num, centr_num, dj, xv, vv);

			if(!status)
			{  
				printf("date error ");
				continue;
			}          
			fprintf(out, "%i\n", pl_num);
			fprintf(out, "%22.12f %22.12f %22.12f\n", xv[0], xv[1], xv[2]);
			fprintf(out, "%22.12f %22.12f %22.12f\n", vv[0], vv[1], vv[2]);
		}
	}

	//printf("Press any key...");
	//getc(stdin);
	fclose(out);
	return 0;
}

int main()
{

	try
	{
		//	CDeEphemeride de("..\\data\\de423.dos");
		//	de.GetCoord(eMERCURY, 2378496.5000000000);
		CheckDEephem();
		//FortEqv();
		//CompareEmpDe();
	}
	catch(const CCustomException &e)
	{
		LPCTSTR stz = e.GetErrorStr();
		int a =0 ;
	}
	catch(const CExcptn &e)
	{
		string stz = e.GetString();
		int a = 0 ;
	}
	return 0;
}
